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We use Rice formulae in order to compute the moments of some level functionals which are 
linked to problems in oceanography and optics: the number of specular points in one and two 
dimensions, the distribution of the normal angle of level curves and the number of dislocations 
in random wavefronts. We compute expectations and, in some cases, also second moments of 
such functionals. Moments of order greater than one are more involved, but one needs them 
whenever one wants to perform statistical inference on some parameters in the model or to test 
the model itself. In some cases, we are able to use these computations to obtain a central limit 
theorem. 
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1. Introduction 

Many problems in applied mathematics require estimations of the number of points, the 
length, the volume and so on, of the level sets of a random function {^(x) :x£ R d }, or 
of some functionals dchned on them. Let us mention some examples which illustrate this 
general situation: 

1 . A hrst example in dimension one is the number of times that a random process 
{X(t) : t £ K} crosses the level u: 

NX(u) = #{seA:X(s)=u}. 

Generally speaking, the probability distribution of the random variable N^(u) is un- 
known, even for simple models of the underlying process. However, there exist some 
formulae to compute E(A^ ) and also higher order moments; see, for example, [6]. 

2. A particular case is the number of specular points of a random curve or a random 
surface. Consider first the case of a random curve. A light source placed at (0, hi) emits a 
ray that is reflected at the point (x, W(x)) of the curve and the reflected ray is registered 
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by an observer placed at (0,h 2 ). Using the equality between the angles of incidence and 
reflection with respect to the normal vector to the curve (i.e., N(x) = (— W'(x),l)), an 
elementary computation gives 

x[r 2 - r\) 



where ctj := hi — W(x) and ri := x 2 + oq, i = 1,2. The points (x,W(x)) of the curve 
such that a; is a solution of (1) are called "specular points". For each Borel subset A of 
the real line, we denote by SP\(A) the number of specular points belonging to A. One 
of our aims is to study the probability distribution of SP\(A). 

3. The following approximation, which turns out to be very accurate in practice for 
ocean waves, was introduced some time ago by Longuct-Higgins ([10, 11]; see also [9]). 
If we suppose that hi and hi are large with respect to W(x) and x, then r% = ai + 
x 2 /(2a.i) + 0(h~ 3 ). (1) can then be approximated by 

, x ai + a 2 x hi + h 2 , , , 1 / 1 1 \ 

W (x) ~ ~ - — — ; — = kx, where k := - h — . (2) 

v ; 2 aia 2 2 hih 2 2\hi h 2 J w 

Set Y(x) := W'(x) — kx and let SP 2 (A) denote the number of roots of Y(x) belonging to 
the set A, an approximation of SPi(A) under this asymptotic. The first part of Section 
2 below will be devoted to obtaining some results on the distribution of the random 
variable SP 2 (R). 

4. Let W : Q C R d -> R d ' with d > d! be a random field and define the level set 

C^(u)={xeQ:W(x) = u}. 

Under certain general conditions, this set is a {d — c?')-dimensional manifold, but, in any 
case, its (d— d')-dimensional Hausdorff measure is well defined. We denote this measure 
by <7d-d'- Our interest will be in computing the mean of the a d-d' -measure of this level 
set, that is, E^-d^Cg^u))], as well as its higher moments. It will also be of interest to 
compute 

Y{s)&a d -d\s) 

» 

where Y(s) is some random field defined on the level set. One can find formulae of this 
type, as well as a certain number of applications, in [5. 14] (d' = 1), [3], Chapter 6, and 
[!]• 

5. Another set of interesting problems is related to phase singularities of random wave- 
fronts. These correspond to lines of darkness in light propagation, or threads of silence in 
sound propagation [4] . In a mathematical framework, they can be defined as the locations 
of points where the amplitudes of waves vanish. If we represent a wave as 



W(x,t)=£(x,t)+iri(x,t), xeK d , 
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where £, 77 are independent homogenous Gaussian random fields, then the dislocations 
are the intersections of the two random surfaces £(x, t) = 0,~/(x, i) = 0. Here, we only 
consider the case d = 2. At fixed time, say t — 0, we will compute the expectation of the 
random variable #{x S S : £(x, 0) = ~/(x, 0) =0}. 

The aim of this paper is threefold: (a) to re-formulate some known results in a modern 
language; (b) to prove a certain number of new results, both for the exact and approxi- 
mate models, especially variance computations in cases in which only first moments have 
been known until now, thus contributing to improve the statistical methods derived from 
the probabilistic results; (c) in some cases, to prove a central limit theorem. 

Rice formulae are our basic tools. For statements and proofs, we refer to the recent book 
[3] . On the other hand, we are not giving full proofs since the required computations are 
quite long and involved; one can find details and some other examples that we do not treat 
here in [2]. For numerical computations, we use MATLAB programs which are available 
at the site http://www.math.univ-toulouse.fr/~azais/prog/programs.html. 

In what follows, A^ denotes the Lebesgue measure in M. d , ad'(B) the d'-dimcnsional 
Hausdorff measure of a Borel set B and M T the transpose of a matrix M. (const) is 
a positive constant whose value may change from one occurrence to another. p%(x) is 
the density of the random variable or vector £ at the point x, whenever it exists. If not 
otherwise stated, all random fields are assumed to be Gaussian and centered. 



2. Specular points in dimension one 

2.1. Expectation of the number of specular points 

We first consider the Longuet-Higgins approximation (2) of the number of SP (x, W(x)), 
that is, 

SP 2 (I) = #{x Gl:Y(x) = W'(x) -kx = 0}. 

We assume that {W(x) :i£R} has C 2 paths and is stationary. The Rice formula for the 
first moment ([3], Theorem 3.2) then applies and gives 



E(SP 2 (J))-= / E(\Y'(x)\\Y(x) = 0) PY{x) (0)dx-- 



\Y'(x)\)-=<p 
v A2 



k.r 



d.r 



G(-fc, v / Al)- 



-,<p 



(3) 



dx, 



where A2 and A4 are the spectral moments of W and 



GO, a) := E{\Z\), Z ~ N{n, a 2 ) = /x[2$(>/er) - 1] + 2<x<p(n/a), 



(4) 



where (p(-) and $(•) are respectively the density and cumulative distribution functions 
of the standard Gaussian distribution. 
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If we look at the total number of specular points over the whole line, we get 



G{k,V\ A ) _ /2A4 1 / ltf T_k* 
71 k\ 2 A4 24 A| 



E( 5P 2 (ffi)) = -IL^ c , A/ _l_ ( 1 + -- + - + , (5) 



which is the result given in [10], part II, formula (2.14), page 846. Note that this quantity 



is an increasing function of jr-. 

We now turn to the computation of the expectation of the number of specular points 
SPi(I) defined by (1). It is equal to the number of zeros of the process {Z(x) := W'(x) — 
mi(x, W(x)) : x £ R}, where 



. x 2 - (hi - w)(h 2 -w) + J[x 2 + (hi - w) 2 }[x 2 + (h 2 - w) 2 ] 

mi(X,w) = — : — : . 

v ' ; x(h l +h 2 -2w) 

Assume that the process { W(x) : x £ R} is Gaussian, centered and stationary, with Ao = 1. 
The process Z is not Gaussian, so we use [3], Theorem 3.4, to get 

E(SPi([a,b]))= dx E(\Z'(x)\\Z(x)=0,W(x)=w) 

J a J —00 

(6) 

V2ti V2ttA 2 
For the conditional expectation in (6), note that 

Z'(x) = W"(x) - ®^( x ,W(x)) - ^{ x ,W{x))W'{x) 

so that under the condition {Z(x) = 0, W(x) ~w}, wc get 

din diiTj 

Z' (x) = W"(x) — K(x,w), where K(x,w) = — — (x,w) H — - — (x,w)mi(x,w). 

ox ow 

Once again, using Gaussian regression, we can write (6) in the form 



E(5P 1 ([a,6])) = ^ V /^t fdx f +X G {mA)e^(-\L 2 + ^^) 



where m = m(x,w) = (A 2 w + K(x,w))/ — A| and G is defined in (4). In (7), the 
integral is convergent as a — > — 00, b — > +00 and this formula is well adapted to numerical 
approximation. 

We have performed some numerical computations to compare the exact expectation 
given by (7) with the approximation (3) in the stationary case. The result depends on 
hi,h%, A4 and A2, and, after scaling, we can assume that A2 = 1. When hi ss h 2 , the ap- 
proximation (3) is very sharp. For example, if h\ — 100, h 2 = 100, A4 = 3, the expectation 
of the total number of specular points over R is 138.2; using the approximation (5), the 
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Figure 1. Intensity of specular points in the case h\ = 100, = 300, A4 = 3. Solid line corre- 
sponds to the exact formula, dashed line corresponds to the approximation (3). 

result with the exact formula is around 2.10 -2 larger (this is the same order as the error 
in the computation of the integral). For h\ = 90, h% = 110, A4 = 3, the results are 136.81 
and 137.7, respectively. If hi = 100, hi = 300, A4 = 3, the results differ significantly and 
Figure 1 displays the densities in the integrand of (6) and (3) as functions of x. 

2.2. Variance of the number of specular points 

We assume that the covariance function E(W(x)W(y)) = T(x — y) has enough regularity 
to perform the computations below, the precise requirements being given in the statement 
of Theorem 1 . 

Writing, for short, S = SP2CR), we have 




(8) 



Using [3], Theorem 3.2, we have 




E{\W"{x) - k\\W"(y) - k\\ W'(x) = kx, W'(y) = ky) 



(9) 



x Pw'(x),W(y)(kx, ky) dxdy, 



where 



Pw(x),w(y)(kx,ky) 
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2tt v / A| -T" 2 (x-y) 



exp 



1 k 2 (X 2 x 2 + 2T" 2 (x - y)xy + X 2 y 2 ) 

2 A|-r' /2 (a;-y) 



(10) 



under the condition that the density (10) does not degenerate for x^-y. 

For the conditional expectation in (9), we perform a Gaussian regression of W"(x) 
(rcsp., W"{y)) on the pair (W'(x), W'(y)). Putting z = x — y, we obtain 



W"(x) 
a y (x) 



9 y (x) + a y (x)W'(x) + b y (x)W'(y), 
T"'(z)T"(z) , , , X 2 T"'{z) 



\l-r» 2 { z y 



by(x) 



\l-r» 2 ( z y 



where y (x) is Gaussian centered, independent of (W'(x), W'(y)). The regression of 
W"{y) is obtained by permuting x and y. 

The conditional expectation in (9) can now be rewritten as an unconditional expecta- 
tion: 



E 



? y (x)-kr'"(z) 



i 



T"{z)x + \ 2 y 



\i-T" 2 (z) 



,(v)-kr»'(-z) 



i 



T"(~z)y + \ 2 x 



(ii) 



Note that the singularity on the diagonal x — y is removable since a Taylor expansion 
shows that for z«0, 



T"'{z) 
It can be checked that 



T"{z)x + \ 2 y 



\l-T»*{z) 



i^(z + 0(z 3 )). 



a 2 {z)=n(o y (x)) 2 )=n{e x {y)) 2 ) = \A 

T"' 2 {z)T"(z) 



\ 2 T" r2 {z) 

\ 2 -r» 2 ( z y 



E(d y (x)9 x (y))=T^(z) 



X 2 -T" 2 (z) 



Moreover, if Xq < +oo, we can show that as z « 0, we have 



o\z). 



1 A2A6 — A| 2 
4 Aa Z 



(12) 

(13) 
(14) 

(15) 



and it follows that the singularity on the diagonal of the integrand in the right-hand side 
of (9) is also removable. 

We will make use of the following auxiliary statement that we state as a lemma for 
further reference. The proof requires some calculations, but is elementary, so we omit it. 
The value of i?(p;0,0) can be found in, for example, [6], pages 211-212. 

Lemma 1. Let 



H(p;n,u)=E(\^ + f i\\r, + u\), 
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where the pair (£,77) is centered Gaussian, E(£ 2 ) = E(?? 2 ) = 1,E(£77) = p. 
Then, if p 2 + v 2 < 1 and < p < 1, 

H{p- p, v) = H{p- 0, 0) + R 2 (p; p, v), 

where 



H(p; 0, 0) = - a/1 - p 2 + — arctan ; = and \R 2 (p; p, v)\ < 3(p 2 + v 2 ). 
n 7T ^1 1 - p 2 

In the next theorem, we compute the equivalent of the variance of the number of spec- 
ular points, under certain hypotheses on the random process W and with the Longuet- 
Higgins asymptotic. This result is new and useful for estimation purposes since it implies 
that, as k — > 0, the coefficient of variation of the random variable S tends to zero at 
a known speed. Moreover, it will also appear in a natural way when normalizing S to 
obtain a central limit theorem. 



Theorem 1. Assume that the centered Gaussian stationary process W = {^(cc) : x £ K} 
is S-dependent, that is, T(z) = if \z\ > 5, and that it has C 4 -paths. Then, as k — > 0, we 
have 



Var(S) = e\ + 0(1), 

rZ 



(16) 



where 



\V2 
1 



2X4 2SX4 



7T V /7T 3 A 2 / 

r"'(z) 2 r"(z) 



j 



+8 9 



a 2 {z)H{p{z)-Q,Q)) 



dz, 



s ^2n{\ 2 + T"{z)) 1 



r< 4 >(*) + 



A 2 -F' 2 (z) 



o~ 2 (z) is defined in (13) and H is defined in Lemma 1. Moreover, as k — > 0, we have 



VVar(5) 
E(S) 



'0k. 



Remarks. 



(1) The 5-dependence hypothesis can be replaced by some weaker mixing condition, 
such as 

|r (i) (z)| < (const)(l + \z\)' a (0<i<4) 
for some a > 1 , in which case the value of 9 should be 



2A 4 1 
7t Jn 



+OO 



a 2 {z)H(p(z);0,0) 1 A 4 
2^\ 2 +r>{z) ny% 



dz. 
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The proof of this extension can be constructed along the same lines as the one we 
give below, with some additional computations. 
(2) The above computations complete the study done in [10] (Theorem 4). In [9], the 
random variable ST^CO is expanded in the Wiener-Hermite chaos. The aforemen- 
tioned expansion yields the same formula for the expectation and also allows a 
formula to be obtained for the variance. However, this expansion is difficult to 
manipulate in order to get the result of Theorem 1. 

Proof of Theorem 1. We use the notation and the computations preceding the state- 
ment of the theorem. 

Divide the integral on the right-hand side of (9) into two parts, corresponding to 
\x — y\ > 5 and \x — y\ < 5, that is, 



E(S(S -1)) 



\x-y\>6 



= h+h- 



\x-y\<6 



(17) 



In the first term, the 5-dependence of the process implies that one can factorize the 
conditional expectation and the density in the integrand. Taking into account that for 
each the random variables W"(x) and W'(x) are independent, we obtain for Ji 



h = 



\x-y\>6 



E(\W"(x) - k\)E(\W"(y) - k\) Pw , {x) (kx) Pw , iy) (ky)dxdy. 



On the other hand, we know that W'{x) (resp., W"(x)) is centered normal with variance 
A2 (resp., A4). Hence, 



1 



\x—y\>6 



2ttA- 



■ exp 



lk 2 (x 2 + y 2 ) 
2 



A, 



da; dy. 



To compute the integral on the right-hand side, note that the integral over the whole x, y 
plane is equal to 1/k 2 so that it suffices to compute the integral over the set \x — y \ < S. 
Changing variables, this last integral is equal to 



+00 



cl.r 



x+S 



1 



-6 



2ttA- 



■ exp 



lk 2 {x 2 + y 2 )1 
2~ 



Xo 



dy 



WA?7t 



0(1), 



where the last term is bounded if k is bounded (remember that we are considering an 
approximation in which k ~ 0). Therefore, we can conclude that 



\x—y\>6 



2ttA- 



■ exp 



lk 2 (x 2 
2" 



' V 



A 2 



da; dy 



k 2 k\/\?n 



0(1), 



from which we deduce, performing a Taylor expansion, that 



h = 



2A 4 



7T 



1 



5 



k 2 kJXoix. 



0(1) 



(18) 
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Let us now turn to I 2 . Using Lemma 1 and the equivalences (12) and (15), whenever 
\z\ = \x — y\ < 5, the integrand on the right-hand side of (9) is bounded by 

{const)[H{p(z);0,0) + k 2 (x 2 +y 2 )]. 

We divide the integral I 2 into two parts. 

First, on the set {(x,y) : \x\ < 26, \x — y\ < 6}, the integral is clearly bounded by some 
constant. 

Second, we consider the integral on the set {(x, y) : x > 26, \x — y \ < 6} . (The symmetric 
case, replacing x > 26 by x < —26, is similar - that is the reason for the factor 2 in what 
follows.) We have (recall that z — x — y) 

7 2 = 0(l) + 2 / / a 2 (z)[H{ P {z);Q,0) + R 2 {p{z);p,v)] 

J J\x-y\<8,x>28 

1 



2n^/\ 2 -T" 2 {z) 

1 k 2 (X 2 x 2 + 2T"(x - y)xy + X 2 y 2 ) 



x exp 



\l-r»(x-y) 



dx dy, 



which can be rewritten as 

7 2 = 0(l) + 2 [ a 2 (z)[H(p(zy,0,0) + R 2 (p(z); t ,,^} 

J-8 



1 



^2n(\ 2 +T"{z)) 

+00 



1 k 2 z 



2,2 



cxp 



I25 yj2n{\ 2 -T"(z)) 
Changing variables, the inner integral becomes 
1 f + °° 1 / 1 



2A 2 -T"(z) VA 2 + r"(z) 2 
2 {x-z/2) 2 



cxp 



-k 



\2-T"(z) 



dx. 



ky/2 



„ ex n-r 2 ) d "i7ii +0(1) 



dz 



(19) 



where t = 2y/2k(26 - z/2)/ v / X 2 - T"(z). 
Substituting this into I 2 , wc obtain 



/ 2 = 0(1) 



kV2' 



To finish, combine (20) with (18), (17), (8) and (5). 



(20) 

□ 
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2.3. Central limit theorem 

Theorem 2. Assume that the process W satisfies the hypotheses of Theorem 1. In ad- 
dition, we assume that the fourth moment of the number of approximate specular points 
on an interval having length equal to 1 is uniformly bounded in k, that is, for all agR 
and < k < 1 , 

E([SP 2 ([a,a + l})f)< (const). (21) 

Then, as k — > 0, 



S - J2\Jnl/k , . , 

=>• A^(0, 1) m distribution. 



Remarks. One can give conditions for the additional hypothesis (21) to hold true. Even 
though they are not nice, they are not costly from the point of view of physical models. 
For example, either one of the following conditions implies (21): 

(i) the paths x ~~> W(x) are of class C 11 (use [3], Theorem 3.6, with m = 4, applied to 
the random process {W'(x) 

(ii) the paths x ~-> W(x) are of class C 9 and the support of the spectral measure has 
an accumulation point (apply [3], Example 3.4, Proposition 5.10 and Theorem 3.4, to 
show that the fourth moment of the number of zeros of W"(x) is bounded). 

Note that the asymptotic here differs from other ones existing in the literature on related 
subjects (compare with, e.g., [7] and [12]). 

Proof of Theorem 2. Let a and f3 be real numbers satisfying the conditions 1/2 < 
a < 1, a + f3 > 1, 2a + /3<2. It suffices to prove the convergence as k takes values on a 
sequence of positive numbers tending to 0. To keep in mind that the parameter is k, we 
use the notation S(k) := S = 5P 2 (K). 

Choose k small enough so that k~ a > 2 and define the sets of disjoint intervals, for 
j = 0, ±1, . . . , ±[fc _/3 ] ([•] denotes integer part), 

= ((J - l)[k- a ]S + 8/2, j [k- a ]5 - S/2), 
I^ = [j[k- a ]S-S/2,j[k- a ]S + 5/2]. 

Each interval Uj has length [fc~"]<5 — 6 and two neighboring intervals Uj are separated 
by an interval of length 6. So, the (^-dependence of the process implies that the random 
variables SP2(Uj),j = 0,±1, . . . , ±[fc~' 3 ], are independent. A similar argument applies to 

SP 2 (l}),j = 0,±l,...,±[k-f>]. 
We write 

T(k)= ]T SP 2 (U*), V k = (Var(S(fc)))" 
\ 3 \<[k-f] 

where the equivalence is due to Theorem 1. 
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The proof is performed in two steps, which easily imply the statement. In the first, it 
is proved that Vk[S(k) — T(k)] tends to in the L 2 of the underlying probability space. 
In the second step, we prove that VkT(k) is asymptotically standard normal. 

Step 1. We first prove that Vk[S(k) — T(k)] tends to in L 1 . Since it is non- negative, 
it suffices to show that its expectation tends to zero. We have 

S(k)-T(k)= Y SP 2 (I^) + Z 1 + Z 2 , 

\j\<[k-0] 

where Z x = SP 2 (-oo, -[Ar /3 ] ■ [k~ a }6 + 8/2), Z 2 = 5P 2 ([/c" /3 ] • [k~ a }6 - 5/2, +oo)). 
Using the fact that E(SP 2 (I)) < (const) Jj tp(kx/^/\ 2 ) dx, we can show that 



V k E(S(k) - T(k)) < (const)k 1 / 2 



^ f£[k- a ]kS\ f + °° 

•'M J J[k- a ][k-P]5 



l£=0 



which tends to zero as a consequence of the choice of a and /?. It suffices to prove that 
V k 2 Var(5(fe) — T(k)) —> as k —> 0. Using independence, we have 

Vax(S(k) - T(k)) = J2 Var(5P 2 (/ J ' £ )) + Var(Z x ) + Var(Z 2 ) 

lil<[fc-' 3 ] 

< nsp 2 (i*)(sp 2 (i*) - 1)) 

\3\<[k~n 

+ E{Z 1 (Z 1 - 1))+E(Z 2 {Z 2 - l))+E(S(k)-T(k)). 

We already know that V k 2 E(S(k) — T(k)) — > 0. Since each Ij can be covered by a fixed 
number of intervals of size one, we know that E(SP 2 (I^)(SP 2 (I^) — 1)) is bounded by a 
constant which does not depend on k and j. Therefore, 

V 2 Y, nSP2(^)(SP 2 ^)-l))<(const)k 1 -^, 
\j\<[k-P] 

which tends to zero because of the choice of /3. The remaining two terms can be bounded 
in a similar form as in the proof of Theorem 1. 

Step 2. T(k) is a sum of independent, but not equidistributed, random variables. To 
prove that it satisfies a central limit theorem, we will use a Lyapunov condition based of 
fourth moments. Set 

Mf := E{[SP 2 {U*) - E(SP 2 (u!f))] m }. 
For the Lyapunov condition, it suffices to verify that 

S~ 4 Y M j^° as A: — s- 0, where E 2 := Y M j- ( 22 ) 

\j\<lk-f>] \j\<lk~f] 
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To prove (22), we divide each interval Uj into p=[k a ] — 1 intervals 1\ , . . . , I p of equal 
size 5. We have 

E(SPi+--- + SP p ) 4 = nSP tl SP l2 SP t3 SP u ), (23) 

l<H,'2,j3,*4<P 

where SPj stands for SP%(Ii) — E(SP2(h))- Since the size of all intervals is equal 
to 5, given the finiteness of fourth moments in the hypothesis, it follows that 
E{SP. tl SPi 2 SPi 3 SP u ) is bounded. 

On the other hand, the number of terms which do not vanish in the sum of the right- 
hand side of (23) is 0(p 2 ). In fact, if one of the indices in (ii, isj *3j *4) differs by more 
than 1 from all the others, then EiSP^SP^SP^SP^) = 0. Hence, 

E[SP 2 (Ulf) - E(SP 2 (Ulf))] 4 < (const)k- 2a 

so that J2\j\<[k-P] Mj = 0(fc~ 2a fc~' 3 ). The inequality 2a + f3 < 2 implies the Lyapunov 
condition. □ 



3. Specular points in two dimensions. 
Longuet-Higgins approximation 

We consider, at fixed time, a random surface depending on two space variables x and y. 
The source of light is placed at (0, 0, hi) and the observer is at (0, 0, /12). The point (x, y) 
is a specular point if the normal vector n{x,y) = {—W x , —W y , 1) to the surface at (x,y) 
satisfies the following two conditions: 

• the angles with the incident ray / = (—x, —y, hi — W) and the reflected ray R = 
{—x,—y,h2 — W) are equal (to simplify notation, the argument (x,y) has been re- 
moved); 

• it belongs to the plane generated by / and R. 

Setting Cii = hi — W and = yj ' x 1 + y 2 + on , i = 1, 2, as in the one-parameter case, we 
have 

u , x a 2 ri-air 2 y a 2 n - a x r 2 , 

W x = a a , W y = . (24) 

x + y A r 2 — n x* +y^ r 2 — r\ 

When h\ and h 2 are large, the system above can be approximated by 

W x = kx, W y = ky, (25) 

under the same conditions as in dimension one. 

Next, we compute the expectation of SP 2 (Q), the number of approximate specular 
points, in the sense of (25), that are in a domain Q. In the remainder of this paragraph, we 
limit our attention to this approximation and to the case in which {VF(x, y) : (x, y) £ R 2 } 
is a centered Gaussian stationary random field. 
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Y(x,y) 



W x (x,y) - kx 
W y (x,y) - ky 



(26) 



Under very general conditions, for example, on the spectral measure of {W(x,y) :x, y £ 
R}, the random field {Y(x,y) :x,y&R} satisfies the conditions of [3], Theorem 6.2, and 
we can write 



E(SP 2 (Q))= / E(\dctY'(x,y)\)pY lx , v) (0)dxdy 



(27) 



since for fixed (x,y), the random matrix Y'(x,y) and the random vector Y(x,y) are 
independent so that the condition in the conditional expectation can be eliminated. The 
density in the right-hand side of (27) has the expression 



Py(x,j/)(0) =P(w m ,w y ){kx,ky) 
1 1 



^ n \J A20A02 — Afi 



cxp 



k 2 



2(A2qAq2 — A n ) 



(A 02 x 2 - 2Xnxy + A 20 y 2 ) 



(28) 



To compute the expectation of the absolute value of the determinant in the right- 
hand side of (27), which does not depend on x,y, we use the method of [4]. Set 
A :=detY'(x,y) = (W xx - k)(W yy - k) - W%. 



We have 



Define 



E(|A|)=E 



1 - cos(At) 
T 2 



dt 



(29) 



h(t) :=E[exp(it[(W M - k)(W yy - k) - WU)]. 



Then 



2 / r +oc 

wo --(J 



l-me[h(t)} 



dt 



(30) 



We now proceed to give a formula for D\e[h(t)}. Define 

/ 1/2 
4=1/2 
\ 0-1 

and denote by S the variance matrix of (W xx , W yy , W x ,y) 



A40 A22 A31 

A22 A04 A13 

. A31 A13 A22 



Rice formulae and Gaussian waves 

Let S 1 / 2 ^!] 1 / 2 = Pdiag(Ai, A 2 , A 3 )P T , where P is orthogonal. Then 

h(t) = e itfe2 E(exp[it((A 1 Z 1 2 - k(s u + s 21 )Z 1 ) + (A 2 Z 2 - k(s 12 + s 22 )Z 2 ) 

+ (A 3 Z 2 -fc( Sl3 + S23 )Z 3 ))]), 
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(31) 



where (Z\, Z 2 , Z 3 ) is standard normal and Sij are the entries of £ 1 / 2 P T . 

One can check that if £ is a standard normal variable and r, [i are real constants, r > 0, 
then 



E(e' 



r(?+M) 2 ) = (1 - 2ir) _1/2 c iT ' i;i/(1_2ir) = 



1 



(l+4 T 2)l/4 



exp 



-2t 



1 + 4t 2 



+ 1 y> + 



l+4r 2 



where ip = h arctan(2r),0 < ip < 7t/4. Substituting this into (31), we obtain 



where, for j = 1, 2, 3: 



dj(t,k) 



^(^■(t)+fc 2 %(i)) 
^=1 > 



(32) 



• dj (t, k) = exp 



k 2 t 2 (gij +s 2 j) 2 
2 l + 4A 2 i 2 



• (i) = - arctan(2Aj t), < ^ < 7t/4; 



1 +2 ( s ii + *aj) a A j 



1 + 4A 2 t 2 



Introducing these expressions into (30) and using (28), we obtain a new formula which 
has the form of a rather complicated integral. However, it is well adapted to numerical 
evaluation. On the other hand, this formula allows us to compute the equivalent as 
k — >• of the expectation of the total number of specular points under the Longuct- 
Higgins approximation. In fact, a first-order expansion of the terms in the integrand 
gives a somewhat more accurate result, one that we now state as a theorem. 



Theorem 3. 



where 



m 2 



E(SP 2 (M 2 )) = ^+0(1), 



(33) 



+ °° 1 ~ [n^i(l + 4A 2 ^ 2 )]-V 2 co S (E 3 3 =i ipjjt)) 
t 2 



dt 



+ x 



1 - 2- 3 / 2 [n, 3 = i(^V T +^7)](l - BiB 2 - B 2 B 3 - B3-B1) 



t 2 



dt 



(34) 
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Figure 2. Intensity function of the specular points for the Jonswap spectrum. 



Aj = Mt) = (1 + 4A 2 jt 2 r 1/2 , B 3 = Bfr) = V /(1-A J )/(1 + ^). 

Note that 7712 depends only on the eigenvalues Ai,A2,A3 and is easily computed 
numerically. We have performed a numerical computation using a standard sea model 
with a Jonswap spectrum and spread function 008(26*). It corresponds to the default 
parameters of the Jonswap function of the toolbox WAFO [13]. The variance matrix of 
the gradient and the matrix S are, respectively, 

1(r *(T 8 °). s ='°- 4 (^ J »)■ 

The integrand in (27) is displayed in Figure 2 as a function of the two space variables 
x,y. The value of the asymptotic parameter 771,2 is 2.52710 -3 . 

We now consider the variance of the total number of specular points in two dimensions, 
looking for analogous results to the one-dimensional case (i.e., Theorem 1), in view of 
their interest for statistical applications. It turns out that the computations become 
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much more complicated. The statements on variance and speed of convergence to zero of 
the coefficient of variation that we give below include only the order of the asymptotic 
behavior in the Longuct-Higgins approximation, but not the constant. However, we still 
consider them to be useful. If one refines the computations, rough bounds can be given on 
the generic constants in Theorem 4 on the basis of additional hypotheses on the random 
field. 

We assume that the real- valued, centered, Gaussian stationary random field {VF(x) : x € 
R 2 } has paths of class C 3 , the distribution of W'(0) does not degenerate (i.e., Var(W(0)) 
is invertible). Moreover, let us consider W"(0), expressed in the reference system xOy 
of R 2 as the 2x2 symmetric centered Gaussian random matrix 



W"(0) 

The function 



W xx {0) W xy {0) 

W X y(0) Wyy(0) 



z ~» A(z) = det[Var(W"(0)z)], 

defined on z = (z\, Z2) T € R 2 , is a non- negative homogeneous polynomial of degree 4 in 
the pair z±, Zi. We will assume the non-degeneracy condition 

min{A(z):||z|| =1} = A>0. (35) 

Theorem 4. Let us assume that {W(x) : x S R 2 } satisfies the above conditions and that 
it is also 8-dependent, 5 > 0, that is, E(W / (x)W (y)) = whenever ||x— y|| > 5. Then, for 
k small enough, 

Var(&P 2 (R 2 )) < ^, (36) 

where L is a positive constant depending on the law of the random field. 

Moreover, for k small enough, by using the result of Theorem 3 and (36), we get 



VVar(ffP 2 (R 2 )) 
E(S/MR 2 )) " 1 ' 

where L\ is a new positive constant. 

Proof. To simplify notation, let us denote T= .ST^R 2 )- We have 

Var(T) = E(T(T - 1)) + E(T) - [E(T)] 2 . (37) 

We have already computed the equivalents as k — > of the second and third term in 
the right-hand side of (37). Our task in what follows is to consider the first term. 

The proof is performed along the same lines as the one of Theorem 1, but instead of 
applying a Rice formula for the second factorial moment of the number of crossings of a 
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one-parameter random process, we need [3], Theorem 6.3, for the factorial moments of a 
2-parameter random field. We have 

E(T(T-1))= / / E(|dctY'(x)||detY'(y)||Y(x)=0,Y(y)=0) 

J Jm 2 xR 2 

x PY(x),Y(y)(0,0)dxdy 



•••dxdy+ / / ••• dxdy = Ji +J 2 . 

x ~y||><5 * ^ll x— yll^^ 

For Ji, we proceed as in the proof of Theorem 1, using the 5-dependence and the 
evaluations leading to the statement of Theorem 3. We obtain 

i - m K 0(1) m) 

One can show that under the hypotheses of the theorem, for small k, one has 

A = ffi>. (39, 

We refer the reader to [2] for the lengthy computations leading to this inequality. In view 
of (37), (33) and (38), this suffices to prove the theorem. □ 



4. The distribution of the normal to the level curve 

Let us consider a modeling of the sea W(x,y,t) as a function of two space variables and 
one time variable. Usual models are centered Gaussian stationary with a particular form 
of the spectral measure /i that is presented, for example, in [3]. We denote the covariance 
byT(x,y,t)=E(W(0,0,0)W(x,y,t)). 

In practice, one is frequently confronted with the following situation: several pictures of 
the sea on time over an interval [0, T] are stocked and some properties or magnitudes are 
observed. If the time T and the number of pictures are large, and if the process is ergodic 
in time, then the frequency of pictures that satisfy a certain property will converge to 
the probability of this property happening at a fixed time. 

Let us illustrate this with the angle of the normal to the level curve at a point "chosen 
at random" . We first consider the number of crossings of a level u by the process W(-,y,t) 
for fixed t and y, defined as 

( u ) = #{ x -°< x < Mi;W{x, y, t) = u}. 

We are interested in computing the total number of crossings per unit time when inte- 
grating over y € [0, M 2 ], that is, 

*7o dt l N ^ t)(u)dy - (40) 
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If the ergodicity assumption in time holds true, then we can conclude that a.s. 

cT pM 2 



1 

T 

where 



dt 



n 



N^f (u) dy - MMN^llT («)) = ^ e " 1/2,l2/A0 °° » 



A a 6c — / A^A y Aj d/x(A a! , A^, At) 

are the spectral moments of VF. Hence, on the basis of the quantity (40), for large T, 
one can make inference about the value of certain parameters of the law of the random 
field. In this example, these are the spectral moments A200 and A oo- 

If two-dimensional level information is available, one can work differently because there 
exists an interesting relationship with Rice formulae for level curves that we explain in 
what follows. We can write (x = (x,y)) 

W'(x, t) = || W(x 3 1) || (cos 8(x, t), sin 6(x, t)) T . 

Using a Rice formula, more precisely, under conditions of [3], Theorem 6.10, 



E 



M-2 



W(-,y,0) 
[0,Mi] 



(u)dy 



= E 


Jc 







|cos6(x,0)|dcri 



/c Q (o,«) 



*2(Q) /A 

Aooo 



71 



200 c _ m 2/(2Aqoo) 



(41) 

where Q = [0,Mi] x [0,M2]. We have a similar formula when we consider sections of the 
set [0, Mi] x [0,M2] in the other direction. In fact, (41) can be generalized to obtain the 
Palm distribution of the angle O. 

Set hg 1 g 2 = I\0. $ 2 ] and, for — 7t < 6\ < 9 2 < tx, define 



F(8 2 ) - F(6»!) := E(<7i({x G Q : W(x, 0) =u;6 x < 6(x, s) < 2 })) 
= E(/ V.e 2 (e(x,s))d ( j 1 (x)ds 



(42) 



cr 2 (Q)E 



|^)((W 2 + (W^ 



cxp(- u 2 /(2A o)) 



OIK) 



Defining A = A200A020 — Ano and assuming o~ 2 (Q) = 1 for case of notation, we readily 
obtain 



F{9 2 )-F{6 X ) 

e -M 2 /(2Aooo) 



(27r)3/2(A)V2yAW 

c -m 2 /(2Aqo) 



K,e2(0)\/^ 2 +y 2 e" (1/(2A))(Ao2 ' T " _2Allxy+A2or) dxdy 



(2tt) 3 / 2 (A+A_) 1 /27a 



nun 
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r+co ,tf 2 / 2 \ 

x / / p exp I — — — — (A+ cos 2 (<^ — k) + A_ sin 2 (95 - k)) \ dp dip, 
Jo Je 1 V 2A + A_ / 

where A_ < A + are the eigenvalues of the covariance matrix of the random vector 
(d x W(Q, 0,0), d y W(0, 0,0)) and k is the angle of the eigenvector associated with 7+. 
Noting that the exponent in the integrand can be written as 1/A_(1 — 7 2 sin 2 (<^ — «)) 
with 7 2 := 1 — A + /A_ and that 



+00 

p 2 exp 



h p 2 \ _ nr 

"V 2H' 



it is easy to obtain that 

F{6 2 ) - F(6 X ) = (const) I (1 - 7 2 sin% - n))'" 1 ' dip. 



From this relation, we get the density g(<p>) of the Palm distribution, simply by dividing 
by the total mass: 

(Wsin 2 ^-*;))- 1 / 2 (l- 7 2 S in 2 (y- K ))-V 2 
■ 9lW /" 7t (l- 7 2 sin 2 (^- K ))-i/2 d ^ 4/C( 7 2 ) ' 1 J 

Here, K. is the complete elliptic integral of the first kind. This density characterizes 
the distribution of the angle of the normal at a point chosen "at random" on the level 
curve. In the case of a random field which is isotropic in (x,y), we have A200 = A020 
and, moreover, Ano = 0, so that g turns out to be the uniform density over the circle 
(Longuet-Higgins says that over the contour, the "distribution" of the angle is uniform 
(cf. [11], page 348)). We have performed the numerical computation of the density (43) 
for an anisotropic process with 7 = 0.5, k = 7t/4. Figure 3 displays the densities of the 
Palm distribution of the angle showing a large departure from the uniform distribution. 



Let us turn to ergodicity. For a given subset Q of R 2 and each t, let us define A 



cr{W(x, y,t):T>t; (x,y) € Q} and consider the cr-algebra of t- invariant events A = f]At- 
We assume that for each pair (x, y), T(x, y, t) — > as t — > +00. It is well known that under 
this condition, the cr-algebra A is trivial, that is, it only contains events having probability 
zero or one (see, e.g., [6], Chapter 7). This has the following important consequence in 
our context. Assume that the set Q has a smooth boundary and, for simplicity, unit 
Lebesgue measure. Let us consider 



Z(t)= / ff(x,t)dCTi(x) (44) 

JC Q (u,t) 

with H(-K,t) = H(W(x,t), VW(x,i)), where WW = (W x , W y ) denotes the gradient in the 
space variables and % is some measurable function such that the integral is well defined. 
This is exactly our case in (42). The process {2(t):tGR} is strictly stationary and, 
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Figure 3. Density of the Palm distribution of the angle of the normal to the level curve in the 
case 7 = 0.5 and k = Tt/4. 

in our case, has a finite mean and is Ricmann-integrable. By the Birkhoff-Khintchinc 
ergodic theorem ([6], page 151), a.s. asT^ +oo, 

where B is the cr-algebra of t- invariant events associated with the process Z(t). Since for 
each t, Z(t) is ^-measurable, it follows that B C A so that E B [Z(0)] = E[Z(0)]. On the 
other hand, the Rice formula yields (taking into account the fact that stationarity of W 
implies that W(0, 0) and \/W(0,0) are independent) 

E[Z(0)] = E[H(u, W(0, 0))|| VW(0, 0)\\}pw{o,a)( u )- 
We consider now the central limit theorem. Let us define 

Z{t) = - [\z(s)-E(Z(0))}ds. (45) 

To compute the variance of 2(t), one can again use the Rice formula for the first moment 
of integrals over level sets, this time applied to the M 2 -valucd random field with parameter 
in K 4 , {(W(xi, s 1 ),W(x 2 ,s 2 )) T ■ (xi,x 2 ) € Q x Q,s 1 ,s 2 € [0,t]} at the level (u,u). We get 

VaxZ(t) = *J (l-j)l(u,s)ds, 
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where 

I(u, s)= [ E[H (xi, 0)H(x 2 , s)\\ VW(xi, 0)|| || W(x 2 , s)|| | W( Xl> 0) = u; VF(x 2 , s) = u] 
Jq 2 

x Pw( Xl ,o),w( X2 M u > u ) dx i dx 2 - {H^u, w (°> 0))||VW(0, 0)||] Pw(Oi0) ( u )) 2 . 

Assuming that the given random field is time-5-dependent, that is, T(x,y,t) = 0\/(x,y) 
whenever t > 8, we readily obtain 



,4 

tV&iZ(t)^2 I(u,s)ds:=cr 2 (u) as i 
Jo 



(46) 



Now, using a variant of the Hocffding-Robbins theorem [8] for sums of (5-dependent 
random variables, we can establish the following theorem. 

Theorem 5. Assume that the random field W and the function H satisfy the conditions 
of [3], Theorem 6.10. Assume, for simplicity, that Q has Lebesgue measure. Then: 

(i) if the covariance j(x,y,t) tends to zero as t—t +oo for every value of (x,y) £ Q, 
we have 

i-t 

Z(«)d«->E[W(«,VW r (0 J 0))||VW(0 ) 0)||]p W (o 1 o)(«). 



where Z(t) is defined by (44)- 

(ii) if the random field W is S-dependent in the sense above, we have 

Vtz(t) ^(o,^ 2 ^)), 

where Z(t) is defined by (45) and o~ 2 {u) by (46). 



5. Application to dislocations of wavefronts 

In this section, we follow the article [4] by Berry and Dennis. Dislocations are lines in space 
or points in the plane where the phase % °f the complex scalar wave ^(x, t) = p(x, t)e Ix ^ x '*^ 
is undefined. With respect to light, they are lines of darkness; with respect to sound, 
threads of silence. Here, we only consider two-dimensional space variables x = (xi,^)- 
It is convenient to express tp by means of its real and imaginary parts: 

V>(x,t)=£(x,t)+i77(x,t). 

Thus, dislocations are the intersection of the surfaces £(x, t) = and i](x,t) = 0. 

Let us quote the authors of [4]: "Interest in optical dislocations has recently revived, 
largely as a result of experiments with laser fields. In low-temperature physics, ip(x.,t) 
could represent the complex order parameter associated with quantum flux lines in a 
superconductor or quantized vortices in a superfluid" (cf. [4] and the references therein) . 
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In what follows, we assume an isotropic Gaussian model. This means that we will 
consider the wavefront as an isotropic Gaussian field 

1 /2 

V;(x, t) = £ exp (i[(k • x) - C |k|t]) (^) dW(k), 

where k = (fci, fc 2 ), |k| = yfcf + fef, II(fc) is the isotropic spectral density and W = (Wi + 
iVF 2 ) is a standard complex orthogonal Gaussian measure on R 2 with unit variance. We 
are only interested in t = and we put £(x) := £(x, 0) and ?7(x) := rj(x, 0). We have, 
setting fc = |k|, 

£(x) = ^ cos((k • x)) V2 d ^(k) - j ^ sin((k • x» (^) V2 dW 3 (k), (47) 

r ? (x) = ^cos((k-x))('ffl) 7 d^ 2 (k)+^sin((k-x))(^^^ 7 dW^k). (48) 

The covariances are 

E[e(x)e(x')]=E[ 7? (x)r ? (x')]= / 9(|x-x'|):= / J (fc|x - x'|)n(fc) dk, (49) 

Jo 

where J v {x) is the Bessel function of the first kind of order v. Moreover, E[£ (ri)?7(r 2 )] = 0. 



5.1. Mean number of dislocation points 

Let us denote by {Z(x) :x€l 2 } a random field having values in R 2 , with coordinates 
£(x),7?(x), which are two independent Gaussian stationary isotropic random fields with 
the same distribution. We are interested in the expectation of the number of dislocation 
points 

d 2 :=E[#{xe5:C(x)=r ? (x)=0}], 

where S is a subset of the parameter space having area equal to 1 . 

Without loss of generality, we may assume that Var(£(x)) = Var(?7(x)) = 1 and for the 
derivatives, we set A 2 = Var(?7j(x)) = Var(£j(x)), i = 1,2. Then, according to the Rice 
formula, 

d 2 =E[|det(Z'(x))|/Z(x) =0]p z( x)(0). 
An easy Gaussian computation gives d 2 = A 2 /(27t) ([4], formula (4.6)). 



5.2. Variance 

Again, let S be a measurable subset of R 2 having Lebesgue measure equal to 1. We have 
Var(iVf (0)) =E(N§ (0)(iVf (0) - 1)) + d 2 - d\ 
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and for the first term, we use the Rice formula for the second factorial moment ([3], 
Theorem 6.3), that is, 

E(iV£(0)(iVf(0)-l)) = / A( Sl ,s 2 )d Sl ds 2 , 
Js 2 

where 

A(s 1 ,s 2 )=E[|dctZ'(s 1 )detZ'(s 2 )||Z(s 1 ) = Z(s 2 ) = 2 ]pz(s 1 ),z(s 2 )(04). 

Here, P denotes the null vector in dimension p. 

Taking into account the fact that the law of the random field Z is invariant under 
translations and orthogonal transformations of R 2 , we have 

A( Sl ,s 2 ) - A((0,0), (r,0)) = A(r) with r = ||si - s 2 ||. 

The Rice function A(r) has two intuitive interpretations. First, it can be viewed as 

A(r) = lim J^E[7V(B((0,0),e)) x N(B((r,0),e))]. 

Second, it is the density of the Palm distribution, a generalization of the horizontal 
window conditioning of the number of zeros of Z per unit surface, locally around the 
point (r, 0), conditionally on the existence of a zero at (0,0) (see [6]). In [4], A(r)jd\ is 
called the "correlation function" . 

To compute A(r), we denote by £i , £ 2 , 771 , r\ 2 the partial derivatives of £,77 with respect 
to first and second coordinate. Therefore, 

A(r) = E[| det Z'(0, 0) det Z'(r, 0) | |Z(0, 0) = Z(r, 0) = 2 ]p ZC o, ) I z(r,o) (0 4 ) 

= E[| (^2 - tamM o)(a% - M(r, 0)||Z(0, 0) = Z(r, 0) = 2 ] (50) 

x Pz(0,0),Z(r,0)(°4)- 

The density is easy to compute: 

1 f°° 
Pz(o,o),z(r,o)(04) = / 27t )2(l - (PM) ' wh ere p(r) = ^ J (fcr)n(fc) dk. 

The conditional expectation turns out to be more difficult to calculate, requiring a long 
computation (we again refer to [2] for the details) . We obtain the following formula (that 
can be easily compared to the formula in [4] since we are using the same notation): 



A(r) 



-4i 



4tt 3 (1-C 2 ) i 2 



1 - 



1 (Z 2 - 2Z 2 i 2 ^ 



(l + t 2 )Z 2 ^/(Z 2 -Z^) 



dt, 
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where we have defined 



C:=p{r), E = p'(r), 
Ai = F { F 



l-C 2 



Z = 



Zi = 



F 2 - H 2 
F 2 

A 2 
1 + Zf 21 



1- F- 



Zn = 



H 

A 2 - 

E 2 C 
l-C 2 

1 + t 2 
1 + Zt 2 ' 



-E/r, F = -p"(r), 

H F(l - C 2 ) - E 2 C 
To F (1~C 2 )-E 2 ' 

2 / E 2 x -a n 

F - 



1-C 2 



Fo 



-p"(o), 
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